setThreadOptions(numThreads=20,stackSize=1e07)

for(modeltype in c("simplified"
                   )){
maindir<-getwd()

source("./auxfiles/loaddata.R")

#a new seed:
set.seed(1514)

setwd(maindir)

testlev<-1


names(upper)<-names(newguess)
names(lower)<-names(newguess)

newguess[which(newguess>upper)]<-upper[which(newguess>upper)]
newguess[which(newguess<lower)]<-lower[which(newguess<lower)]

#Shutting down PT work moments:
newguess[grepl("part",names(newguess))]<-0
upper[grepl("part",names(newguess))]<-0
lower[grepl("part",names(newguess))]<-0

origuess<-newguess
paramnums<-which(upper!=lower)


cfacts<-NULL
#Size of policy adjustment for counterfactuals
marginval <- 1.1


particles<-mpi.universe.size()-1 

param_var<-as.matrix(read.csv(paste0("./modeloutput/newguess_varcov_",modeltype,".csv"),header=FALSE))
paramdraws<-list()
for(i in 1:(particles)){
paramdraws[[i]]<-rmvnorm(100,origuess[paramnums],sigma=param_var)
for(j in 1:nrow(paramdraws[[i]])){
  paramdraws[[i]][j,c("wagevar","spwagevar")]<-pmax(paramdraws[[i]][j,c("wagevar","spwagevar")],1e-4)
  paramdraws[[i]][j,"wagecorr"] <- pmax(paramdraws[[i]][j,"wagecorr"],0)
  paramdraws[[i]][j,"wagecorr"] <- pmin(paramdraws[[i]][j,"wagecorr"],0.99)
  paramdraws[[i]][j,grepl("pi1.|pi2.|pi0.",colnames(paramdraws[[i]]))] <- pmax(paramdraws[[i]][j,grepl("pi1.|pi2.|pi0.",colnames(paramdraws[[i]]))],0)
  paramdraws[[i]][j,grepl("pi1.|pi2.|pi0.",colnames(paramdraws[[i]]))] <- pmin(paramdraws[[i]][j,grepl("pi1.|pi2.|pi0.",colnames(paramdraws[[i]]))],1)
}
}

shocks.wagenoiseraw <- shocks.wagenoise/wagepars$noisevar
shocks.spwagenoiseraw <- shocks.spwagenoise/spwagepars$noisevar

setThreadOptions(numThreads=10,stackSize=1e07)

#########################
#PARALLEL STUFF FOR LINUX:
#########################
print(mpi.universe.size())
#numThreads<-10
seed=6327
partic_mult<-1
##partic_mult<-4
##threads<-10
#########################


cl <- makeCluster(particles, type="MPI",outfile="")
registerDoSNOW(cl)
clusterEvalQ(cl,{
  library(Rcpp)
  library(data.table)
  library(RcppArmadillo)
  library(RcppParallel)
  library(ggplot2)
  library(parallel)
  library(snow)
  library(doSNOW)
  library(Rmpi)
  library(spousalVFIParallel)
  library(unixtools)
  library(readr)
  library(fixest)
  library(readstata13)
  
  source("./auxfiles/plotrules.R")
  source("./auxfiles/evalfunctions.R")
  source("Welfare_ParallelSpouseC.R")}
)
#########################

ex<- Filter(function(x) is.function(get(x,.GlobalEnv)),ls(.GlobalEnv))
clusterExport(cl,ex)

clusterEvalQ(cl,{
  setThreadOptions(numThreads=10,stackSize=1e07)
}
)
#########################


#foreach(particle = 1:(particles))%dopar%{
  foreach(particle = 1:(particles+1))%dopar%{
  Cfactmat_ss <- NULL
  Cfactmat_fs <- NULL
  Cfactmat_stringency <- NULL
  Cfactmat_marstringency <- NULL
  Cfactmat_singstringency <- NULL
  Cfactmat_earlydi <- NULL
  Cfactmat_meanstest <- NULL
  Cfactmat_meanstest20k <- NULL
  Cfactmat_meanstest40k <- NULL
  Cfactmat_meanstest90k <- NULL
  Cfactmat_ss_meanstest <- NULL
  Cfactmat_fs_meanstest <- NULL
  Cfactmat_stringency_meanstest <- NULL
  Cfactmat_marstringency_meanstest <- NULL
  Cfactmat_singstringency_meanstest <- NULL
  Cfactmat_earlydi_meanstest <- NULL
  Cfactmat_fastdi <- NULL
  Cfactmat_fastdi_meanstest <- NULL
  
  mydraws<-paramdraws[[particle]]
  for(n in 1:nrow(paramdraws[[particle]])){
  newguess<-origuess
  newguess[paramnums]<-mydraws[n,]
    
  ###### Upper bound #########
  temp<-wagesetup_fit(newguess,0,0,wage_gridsize=10,spwage_gridsize=5#,filesuffix=typename
  )
  #Wage<-temp$Wage
  #spWage<-temp$spWage
  wagepars<-temp$wagepars
  spwagepars<-temp$spwagepars
  rm(temp)
  
  shocks.wagenoise <- shocks.wagenoiseraw*wagepars$noisevar
  shocks.spwagenoise <- shocks.spwagenoiseraw*wagepars$spnoisevar
  
  shocks.wage <- exp(shocks.wageraw*wagepars$wagevar)
  shocks.spwage <- exp((shocks.spwageraw*(1-wagepars$wagecorr^2)^0.5 + shocks.wageraw*wagepars$wagecorr)*spwagepars$wagevar)
  shocks.wagereal <- exp(shocks.wageraw*wagepars$wagevar)
  shocks.spwagereal <- exp((shocks.spwageraw*(1-wagepars$wagecorr^2)^0.5 + shocks.wageraw*wagepars$wagecorr)*spwagepars$wagevar)
  
  temp<-paramsetup(newguess,wagepars,spwagepars,
             ptwork_halfcost=ptwork_halfcost,
             single_samecost=single_samecost,
             ptwork_halfdis=ptwork_halfdis,
             single_samedis=single_samedis,
             spouse_samejobrisk=spouse_samejobrisk,
             single_sameprefs = single_sameprefs,
             health_simpleprefs = health_simpleprefs,
             single_sameallowanceprob=single_sameallowanceprob,
             work_noallowanceprob=work_noallowanceprob)
  params<-temp$params
  param_input<-temp$param_input
  rm(temp)
  
  
    simdata<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                        shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                        marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  simdata<-cleansim(simdata,select=select,selectC=selectC)
print(paste0("base done for particle ",particle," draw ",n))

uflow<-((simdata[age<62,mean(flow)]*(1-params$gamma))^(1/(1-params$gamma)))

#First, distribution of compensating variations. This one is easy, since it is a function of parameters:
  CVs<-c(cost(newguess[c("theta_sev")]/2,uflow),
    cost(newguess[c("sptheta")],uflow),
    cost(newguess[c("eta")],uflow),
    cost(newguess[c("speta")],uflow),
    cost(newguess[c("eta_sev")]/2,uflow),
    cost(newguess[c("speta_mod")],uflow),
    param_input["F3"]/2
    )
  CVmat <- rbind(NULL,c(CVs,particle))
  write_csv(x=as.data.frame(CVmat),path=paste0("./modeloutput/bootstrap/bootstrap_CV_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  

  skip<-1
  if(skip==0){
    
  ##SS generosity
  params_ss<-params
  params_ss$SSgenerosity <- marginval
  simdata_ss<-Valsim_solve(as.vector(param_input),params_ss,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                  shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                  marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  simdata_ss<-cleansim(simdata_ss,select=select,selectC=selectC)
  print(paste0("ss done for particle ",particle," draw ",n))
  cfacts_ss1 <- cfacts(simdata, simdata_ss, params)
  print(paste0("ss fully done for particle ",particle," draw ",n))
  
  #FS generosity
  params_fs<-params
  params_fs$FSgenerosity <- marginval
  simdata_fs<-Valsim_solve(as.vector(param_input),params_fs,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                           shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                           marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  print(paste0("fs done for particle ",particle," draw ",n))
  simdata_fs<-cleansim(simdata_fs,select=select,selectC=selectC)
  cfacts_fs1 <- cfacts(simdata, simdata_fs, params)
  print(paste0("fs fully done for particle ",particle," draw ",n))
  
  #Stringency
  param_input_stringency<-param_input
  param_input_stringency[grepl("pi1.|pi2.|pi0.",names(param_input))]<-pmin(param_input_stringency[grepl("pi1.|pi2.|pi0.",names(param_input))]*marginval,1)
  simdata_stringency<-Valsim_solve(as.vector(param_input_stringency),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                           shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                           marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  print(paste0("string done for particle ",particle," draw ",n))
  simdata_stringency<-cleansim(simdata_stringency,select=select,selectC=selectC)
  cfacts_stringency1 <- cfacts(simdata, simdata_stringency, params)
  print(paste0("string fully done for particle ",particle," draw ",n))
  
  #Stringency, singles only
  param_input_singstringency<-param_input
  param_input_singstringency[grepl("pi1.|pi2.|pi0.",names(param_input)) & grepl("single",names(param_input))]<-pmin(param_input_singstringency[grepl("pi1.|pi2.|pi0.",names(param_input))  & grepl("single",names(param_input))]*marginval,1)
  simdata_singstringency<-Valsim_solve(as.vector(param_input_singstringency),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                       shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                       marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  simdata_singstringency<-cleansim(simdata_singstringency,select=select,selectC=selectC)
  cfacts_singstringency1 <- cfacts(simdata, simdata_singstringency, params)
  
  #Stringency, married only
  param_input_marstringency<-param_input
  param_input_marstringency[grepl("pi1.|pi2.|pi0.",names(param_input)) & grepl("married",names(param_input))]<-pmin(param_input_marstringency[grepl("pi1.|pi2.|pi0.",names(param_input))  & grepl("married",names(param_input))]*marginval,1)
  simdata_marstringency<-Valsim_solve(as.vector(param_input_marstringency),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                      shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                      marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT)
  simdata_marstringency<-cleansim(simdata_marstringency,select=select,selectC=selectC)
  cfacts_marstringency1 <- cfacts(simdata, simdata_marstringency, params)
  
  
  #Early DI
  simdata_earlydi<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                           shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                           marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT, earlyDI=2)
  simdata_earlydi<-cleansim(simdata_earlydi,select=select,selectC=selectC)
  cfacts_earlydi1 <- cfacts(simdata, simdata_earlydi, params)
  
  #Means test - 20k
  simdata_meanstest<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                  shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                  marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT,DI_meanstest=20000)
  simdata_meanstest<-cleansim(simdata_meanstest,select=select,selectC=selectC)
  cfacts_meanstest20k1 <- cfacts(simdata, simdata_meanstest, params)

  #Means test - 40k
  simdata_meanstest<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                  shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                  marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT,DI_meanstest=40000)
  simdata_meanstest<-cleansim(simdata_meanstest,select=select,selectC=selectC)
  cfacts_meanstest40k1 <- cfacts(simdata, simdata_meanstest, params)

  #Means test - 90k
  simdata_meanstest<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                  shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                  marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT,DI_meanstest=90000)
  simdata_meanstest<-cleansim(simdata_meanstest,select=select,selectC=selectC)
  cfacts_meanstest90k1 <- cfacts(simdata, simdata_meanstest, params)
  
  
  #Means test (main)
  simdata_meanstest<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                  shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                  marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT,DI_meanstest=testlev)
  simdata_meanstest<-cleansim(simdata_meanstest,select=select,selectC=selectC)
  cfacts_meanstest1 <- cfacts(simdata, simdata_meanstest, params)
  
  ##SS generosity and means test
  params_ss<-params
  params_ss$SSgenerosity <- marginval
  simdata_ss_meanstest<-Valsim_solve(as.vector(param_input),params_ss,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                     shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                     marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT, DI_meanstest=testlev)
  simdata_ss_meanstest<-cleansim(simdata_ss_meanstest,select=select,selectC=selectC)
  cfacts_ss_meanstest1 <- cfacts(simdata_meanstest, simdata_ss_meanstest, params)
  
  #Stringency  and means test
  param_input_stringency<-param_input
  param_input_stringency[grepl("pi1.|pi2.|pi0.",names(param_input))]<-pmin(param_input_stringency[grepl("pi1.|pi2.|pi0.",names(param_input))]*marginval,1)
  simdata_stringency_meanstest<-Valsim_solve(as.vector(param_input_stringency),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                             shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                             marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT, DI_meanstest=testlev)
  simdata_stringency_meanstest<-cleansim(simdata_stringency_meanstest,select=select,selectC=selectC)
  cfacts_stringency_meanstest1 <- cfacts(simdata_meanstest, simdata_stringency_meanstest, params)
  
  #Stringency (single only)  and means test
  param_input_singstringency<-param_input
  param_input_singstringency[grepl("pi1.|pi2.|pi0.",names(param_input)) & grepl("single",names(param_input))]<-pmin(param_input_singstringency[grepl("pi1.|pi2.|pi0.",names(param_input))  & grepl("single",names(param_input))]*marginval,1)
  simdata_singstringency_meanstest<-Valsim_solve(as.vector(param_input_singstringency),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                                 shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                                 marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT, DI_meanstest=testlev)
  simdata_singstringency_meanstest<-cleansim(simdata_singstringency_meanstest,select=select,selectC=selectC)
  cfacts_singstringency_meanstest1 <- cfacts(simdata_meanstest, simdata_singstringency_meanstest, params)
  
  #FS generosity  and means test
  params_fs<-params
  params_fs$FSgenerosity <- marginval
  simdata_fs_meanstest<-Valsim_solve(as.vector(param_input),params_fs,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                     shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                     marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT, DI_meanstest=testlev)
  simdata_fs_meanstest<-cleansim(simdata_fs_meanstest,select=select,selectC=selectC)
  cfacts_fs_meanstest1 <- cfacts(simdata_meanstest, simdata_fs_meanstest, params)
  
  
  #Stringency (married only)  and means test
  param_input_marstringency<-param_input
  param_input_marstringency[grepl("pi1.|pi2.|pi0.",names(param_input)) & grepl("married",names(param_input))]<-pmin(param_input_marstringency[grepl("pi1.|pi2.|pi0.",names(param_input))  & grepl("married",names(param_input))]*marginval,1)
  simdata_marstringency_meanstest<-Valsim_solve(as.vector(param_input_marstringency),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                                shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                                marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT, DI_meanstest=testlev)
  simdata_marstringency_meanstest<-cleansim(simdata_marstringency_meanstest,select=select,selectC=selectC)
  cfacts_marstringency_meanstest1 <- cfacts(simdata_meanstest, simdata_marstringency_meanstest, params)
  
  
  #Early DI  and means test
  simdata_earlydi_meanstest<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
                                          shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
                                          marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT, earlyDI=2, DI_meanstest=testlev)
  simdata_earlydi_meanstest<-cleansim(simdata_earlydi_meanstest,select=select,selectC=selectC)
  cfacts_earlydi_meanstest1 <- cfacts(simdata_meanstest, simdata_earlydi_meanstest, params)
  
  # #Fast DI
  # simdata_fastdi<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
  #                               shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
  #                               marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT, earlyDI=1)
  # simdata_fastdi<-cleansim(simdata_fastdi,select=select,selectC=selectC)
  # cfacts_fastdi1 <- cfacts(simdata, simdata_fastdi, params)
  # 
  # #Fast DI  and means test
  # simdata_fastdi_meanstest<-Valsim_solve(as.vector(param_input),params,healthprob, marriageprob,foodstamps,numkids,oecd,aperm(Asset),aperm(Wage), aperm(spWage),
  #                               shocks.wage,shocks.spwage,shocks.health,shocks.allow,shocks.reassess,shocks.jobloss,shocks.spjobloss, shocks.marriage,
  #                               marriageprob_init,numsims,lumpsum=0,print=0,workDI=0,onlyFT=onlyFT, earlyDI=1, DI_meanstest=testlev)
  # simdata_fastdi_meanstest<-cleansim(simdata_fastdi_meanstest,select=select,selectC=selectC)
  # cfacts_fastdi_meanstest1 <- cfacts(simdata, simdata_fastdi_meanstest, params)


  
  Cfactmat_ss <- rbind(NULL,c(cfacts_ss1,particle))
  Cfactmat_fs <- rbind(NULL,c(cfacts_fs1,particle))
  Cfactmat_stringency <- rbind(NULL,c(cfacts_stringency1,particle))
  Cfactmat_singstringency <- rbind(NULL,c(cfacts_singstringency1,particle))
  Cfactmat_marstringency <- rbind(NULL,c(cfacts_marstringency1,particle))
  Cfactmat_earlydi <- rbind(NULL,c(cfacts_earlydi1,particle))
  Cfactmat_meanstest20k <- rbind(NULL,c(cfacts_meanstest20k1,particle))
  Cfactmat_meanstest40k <- rbind(NULL,c(cfacts_meanstest40k1,particle))
  Cfactmat_meanstest90k <- rbind(NULL,c(cfacts_meanstest90k1,particle))
  Cfactmat_meanstest <- rbind(NULL,c(cfacts_meanstest1,particle))

  Cfactmat_ss_meanstest <- rbind(NULL,c(cfacts_ss_meanstest1,particle))
  Cfactmat_fs_meanstest <- rbind(NULL,c(cfacts_fs_meanstest1,particle))
  Cfactmat_stringency_meanstest <- rbind(NULL,c(cfacts_stringency_meanstest1,particle))
  Cfactmat_singstringency_meanstest <- rbind(NULL,c(cfacts_singstringency_meanstest1,particle))
  Cfactmat_marstringency_meanstest <- rbind(NULL,c(cfacts_marstringency_meanstest1,particle))
  Cfactmat_earlydi_meanstest <- rbind(NULL,c(cfacts_earlydi_meanstest1,particle))
  
  write_csv(x=as.data.frame(Cfactmat_ss),path=paste0("./modeloutput/bootstrap/bootstrap_ss_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_fs),path=paste0("./modeloutput/bootstrap/bootstrap_fs_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_stringency),path=paste0("./modeloutput/bootstrap/bootstrap_stringency_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_singstringency),path=paste0("./modeloutput/bootstrap/bootstrap_singstringency_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_marstringency),path=paste0("./modeloutput/bootstrap/bootstrap_marstringency_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_earlydi),path=paste0("./modeloutput/bootstrap/bootstrap_earlydi_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_meanstest20k),path=paste0("./modeloutput/bootstrap/bootstrap_meanstest20k_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_meanstest40k),path=paste0("./modeloutput/bootstrap/bootstrap_meanstest40k_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_meanstest90k),path=paste0("./modeloutput/bootstrap/bootstrap_meanstest90k_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_meanstest),path=paste0("./modeloutput/bootstrap/bootstrap_meanstest_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  
  write_csv(x=as.data.frame(Cfactmat_ss_meanstest),path=paste0("./modeloutput/bootstrap//bootstrap_ss_meanstest_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_fs_meanstest),path=paste0("./modeloutput/bootstrap/bootstrap_fs_meanstest_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_stringency_meanstest),path=paste0("./modeloutput/bootstrap/bootstrap_stringency_meanstest_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_singstringency_meanstest),path=paste0("./modeloutput/bootstrap/bootstrap_singstringency_meanstest_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_marstringency_meanstest),path=paste0("./modeloutput/bootstrap/bootstrap_marstringency_meanstest_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  write_csv(x=as.data.frame(Cfactmat_earlydi_meanstest),path=paste0("./modeloutput/bootstrap/bootstrap_earlydi_meanstest_",modeltype,".csv"),append=TRUE,col_names=FALSE)
  
  }

#  Cfactmat_fastdi_meanstest <- cbind(Cfactmat_fastdi_meanstest,(cfacts_fastdi_meanstest1-cfacts_fastdi_meanstest0)/(param1-param0))
#  Cfactmat_fastdi <- cbind(Cfactmat_fastdi,(cfacts_fastdi1-cfacts_fastdi0)/(param1-param0))
#  write_csv(x=as.data.frame(Cfactmat_fastdi),path=paste0("gradient_fastdi_",modeltype,".csv"),append=FALSE,col_names=FALSE)
#  write_csv(x=as.data.frame(Cfactmat_fastdi_meanstest),path=paste0("gradient_fastdi_meanstest_",modeltype,".csv"),append=FALSE,col_names=FALSE)
  }
}
}



